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ABSTRACT 

We develop a formulation for constructing and examining rapidly rotating Newtonian 
neutron star models that contain two superfluids, taking account of the effect of the 
' rotation velocity difference between two superfluids. We assume neutron stars to be 

, composed of the superfluid neutrons and the mixture of the superfluid protons and 

' the normal fluid electrons. To describe Newtonian dynamics of the two superfluids, 

the Newtonian version of the so-called two-fluid formalism is employed. The effect 
of the rotation velocity difference on the structure of equilibrium state is treated as 
, a small perturbation to rapidly rotating superfluid stars whose angular velocities of 

two superfluids are assumed to be exactly the same. We derive basic equations for 
the perturbed structures of rapidly rotating superfluid stars due to the rotation ve- 
locity difference between two superfluids. Assuming the superfluids to obey a simple 
Q , analytical equation of state proposed by Prix, Comer, and Andersson, we obtain nu- 

' merical solutions for the perturbations and find that the density distributions of the 

, superfluids are strongly dependent on the parameter a which appears in the analyti- 

cal equation of state and characterizes the so-called symmetry energy, ft is also found 
that if Prix et al.'s analytical equation of state is assumed, the perturbations can be 
' represented in terms of the universal functions that are independent of the parameters 

^ \ of the equation of state. 
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1 INTRODUCTION 

To investigate properties of equilibrium configurations of rotating neutron stars, so far, most neutron star models have been 
obtained by assuming neutron star matter to be a one-constituent perfect fluid (for a review, see, e.g. Stergioulas 2003) . This 
treatment of equilibrium states of neutron stars seems to be quite reasonable as a first approximation for examining global 
properties of neutron stars such as the gravitational mass, the radius, or the maximum rotation frequency. However, it has long 
been suggested that neutrons in the inner crust and neutrons and protons in the core of neutron stars are in superfluid states 
when the interior temperatures cool down below Tc ~ lO^K (e.g., Shapiro & Teukolsky 1983). Since the interior temperature 
of neutron stars is believed to cool down quickly via the neutrino emission (e.g., Baym & Pethick 1979), it is likely that 
each of many observable neutron stars, except newly born ones, has a core containing superfluids. Although superfluidity 
in the interior might be one of the important ingredients that affect the structures of neutron stars, the superfluidity has 
been neglected in most studies on equilibrium configurations of rotating neutron stars. Thus, it is necessary to examine how 
superfluidity affects fundamental properties of rotating neutron stars. 

The superfluidity in neutron stars has been mainly argued in connection with the gfitch phenomenon and the post- 
glitch relaxation observed in pulsars. The glitch is a sudden decrease of an observed pulse period Tp of a pulsar, whose largest 
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magnitude is ATp/Tp ^ —10"'', that is followed by a period of continuous relaxation. In most explanations of glitch phenomena, 
a two-component model has been employed to describe rotational dynamics of neutron stars containing superfluid neutrons. 
In those treatments, one of the components is a mixture of the charged particles such as protons (or ions) and electrons, which 
are assumed to co-rotate because the charged components are strongly coupled each other due to short-range electromagnetic 
interactions (Easson 1979, Alpar et al. 1984), and the other the superfluid neutrons, which are supposed to rotate separately 
from the charged components because the superfluid neutrons arc inviscous and loosely couple to the charged components. In 
this paper, we call the mixture of the charged components "protons" for brevity. It has been shown that those two-component 
superfluid models have succeeded in explaining observed features of the glitches and the post-glitch relaxations qualitatively 
(Anderson & Itoh 1975; Alpar et al. 1984; Sedrakian et al. 1995, Link & Epstein 1996). Those results mean that the existence of 
superfluid components inside neutron stars is supported not only from the theory of nuclear physics but also from observations 
of pulsars. 

If they co-rotate, as mentioned in Prix & Rieutord (2002), two constituents in a rotating superfluid star behave as one- 
constituent ordinary fluid, because the so-called entrainment effect, due to which the momentum of one superfluid constituent 
is dependent on the mass current of all superfluid constituents, does not operate and two constituents are in a chemical 
equilibrium state. On the other hand, specific characteristics of neutron stars due to the superfluidity are considered to be 
observed when two constituents rotate at different rotation rates. To study the effect of rotation velocity difference between 
two constituents on the global structures of rotating neutron stars with superfluidity, Prix (1999) derived basic equations for 
the equilibrium conflgurations and obtained analytical solutions for slowly rotating superfluid stars within the framework of 
Newtonian dynamics, generalizing the so-called Chandrasekhar-Milne expansion for ordinary fluid stars (see, e.g. Tassoul 1978), 
in which rotational effects are treated as perturbations to non-rotating spherical stars. To describe dynamics of superfluids, 
he made use of a variant of the so-called two-fluid formalism, which had been mainly developed by Carter, Langlois, and 
their co-workers (Carter 1989; Carter & Langlois 1998; Langlois, Sedrakian, & Carter 1998; Prix 2002). In order to obtain 
analytical solutions, Prix (1999) employed an equation of state whose functional form is P oc p^, where P and p axe the 
pressure and density, respectively. More recently, Prix, Comer, & Andersson (2002a) extended the Prix formalism so as to 
include the entrainment effect between two superfluids in order to investigate its effect on the stellar structures. In the same 
paper, Prix et al. (2002a) also argued the effect of the so-called symmetry energy, which can be included as a parameter in the 
equation of state, and found that the effect of symmetry energy is significant in determination of the density distributions of 
the neutrons and the protons. As for relativistic rotating superfluid stars, Andersson & Comer (2001) calculated equilibrium 
configurations of slowly rotating superfluid neutron stars, extending the slow rotation approximation formalism devised by 
Hartle (1967) (sec, also Hartlo & Thorne 1968) for ordinary fluid stars to relativistic two-constituent superfluid stars. Very 
recently, Prix, Novak, & Comer (2002b) have obtained preliminary results for rapidly rotating superfluid stars, without the 
slow rotation approximation. 

The purpose of this paper is to improve our understanding of properties of equilibrium configurations of rapidly rotating 
neutron stars with superfiuids. We calculate equilibrium configurations of rapidly rotating neutron stars with neutron and 
proton superfluids. In this paper, we are concerned with rotating neutron stars of two superfluids whose angular velocities are 
different. We however assume that the rotation velocity difference between two superfluids is very small in comparison with the 
angular velocity of the star. This assumption could be reasonable for models of superfluid neutron stars because observations 
of pulsar glitches show that the rotation velocity of the protons may differ from that of the neutrons, but the amount of the 
rotation velocity difference is not very large. We can therefore treat the effect of the rotation velocity difference between two 
superfluids on the structures of rapidly rotating neutron stars as a small perturbation to equilibrium configurations of rapidly 
rotating neutron stars in which two superfluids co-rotate. For the dynamics of superfluids inside stars, we employ a Newtonian 
version of the two-constituent formalism developed by Prix (2002). In order to obtain basic equations for determining the 
effect of the rotation velocity difference between two superfluids, equations of superfluid hydrostatic equilibrium and the 
Poisson equation are expanded in terms of the rotation velocity difference between two superfluids. To obtain the numerical 
solutions of equilibrium configurations, we make use of a variant of the so-called Self Consistent Field method (SCF) devised 
by Ostriker & Mark (1968). The equation of state for the superfluid we use is the analytical one employed by Prix et al. 
(2002a), which is a natural extension of an A'" = 1 polytropic equation of state for a barotropic ordinary fluid to the case of two 
superfluids. In §2 we present the basic equations employed in this paper for the rapidly rotating superfluid stars. Numerical 
results are given in §3, and §4 is devoted for summary and discussion. 

2 FORMULATION 

2.1 Two-constituent formalism for Newtonian superfluid dynamics 

In this paper, a neutron star is assumed to be composed of superfluid neutrons and a mixture of superfluid protons and 
normal fluid electrons, because the interior temperatures of old neutron stars are much lower than the transition temperature 
to neutron and proton superfluids, which is considered to be Tc ~ 10^ K (see, e.g., Epstein 1988). We further assume that 
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this mixture of charged particles is perfectly in the state of charge neutrality because of strong coupling between the protons 
and the electrons by electromagnetic interactions. The electrons therefore co-move with the protons, and their motions can be 
described with one fluid velocity Vp. Through this paper, we call the mixture of the the protons and the electrons "protons" 
for brevity as mentioned before. To describe the dynamics of the superfluids, we make use of the two-constituent formalism 
for Newtonian superfluid dynamics, which has been recently developed by Prix (2002). 

We take account of no "transfusion" between the neutrons and the protons. The particle numbers of the neutron and the 
proton must be therefore conserved separately. The particle number conservation equations for the neutron and the proton 
are then given by 

Va{n„v^) = 0, \/ainpV;)=0, (1) 

where n„ and rip are the number densities of the neutron and the proton, and uJJ and Vp denote the fluid velocities of the 
neutron and the proton, respectively. Here, Va means the covariant derivative in the three dimensional flat space. The mass 
densities of two superfluids can be written as 

p„ — mn„ , pp — mup , (2) 

where m is the mass of the neutron, and we assume the mass of the proton to be equal to that of the neutron. The fundamental 
quantity of the two-constituent formalism for the Newtonian superfluid dynamics is the internal energy density £, which is 
a function of p„, pp, and A^, where — A^A^ and A" — Vp — v"^. Note that we have neglected the entropy s carried 
by the normal fluid of the electron for simplicity. This treatment is justifled for neutron stars whose internal temperatures 
are sufficiently low. This function £ defines several basic thermodynamical quantities which describe the dynamics of two 
superfluids in terms of the total difference as follows 

d£ = jl„ dpn + p.p dpp + a dA^ , (3) 

where jln and fip denote the specific chemical potentials of the neutron and the proton, respectively. The function a represents 
the strength of the entrainment effect. By using a variant of the Gibbs-Duhem relation, the generalized pressure of two fluids 
is defined by 

dP = p„ dfin + Pp dpp — a dA^ , (4) 
and we can obtain the relation, 

P = Pn fj.n + Pp fj-p — £ ■ (5) 

According to Prix (2002), the Euler equations for two superfluids can be given by 

{dt + C^Jiv„,a + En A,) + Va (^An + $ - i \v„f^ = , (6) 
{dt + C^^){Vp,a - EpAa) + Va (/ip + $ - i = , (7) 

where En = 2a/ pn and Sp = 2a/ pp. Here, Cu means the Lie derivative along the vector field m". Note that we have considered 
no direct interaction force between two superfiuids such as the mutual friction or the Magnus-type force. In equations ® and 
Q, the function $ is the gravitational potential, which is determined with the total mass density p, given by p = p„ + pp, 
through the Poisson equation 

V„V"$ = 47rGp, (8) 

where G is the gravitational constant. Considering the exterior differentiation of equations ® and ||7|l, we can obtain the 
vorticity equations, given by 

(l9t + 'C„„)LJ„,ab = , {dt + C^^)uJp^ab = , (9) 

where uin,ab and uop^ab are, respectively, the exterior derivatives of the one-forms v^^a + enA^ and Vp^a — £pAa, which are 
defined by 

t^n,ab = 9a(Un,b + EnAft) - dbiVn,a + EriAa) , UJp^ab = da{Vp,b — EpAf,) - db{Vp,a " EpAa) . (10) 



2.2 Stationary and axisymmetric equilibrium configurations 

We consider stationary and axisymmetric equilibrium configurations. Thus, the time and azimuthal derivatives of physical 
quantities must vanish because a star is assumed to be in a stationary and axisymmetric state. In this paper, we further 
assume that the neutrons and the protons, respectively, rotate with the angular velocities n„ and Q,p around the axis of 
rotation, and that no meridional circulation is present in a superfluid star. The fluid velocities for two fluids are then given by 
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= ^n^" , Vp = ripVp" , (11) 

where ip"" is the rotational Kilhng vector. In this paper, we employ the spherical polar coordinates (r,6,ip). The components 
of (/j" can be written as y?" = in the spherical polar coordinates. First of all, let us consider the vorticity conservation for 
the flow given by equation lll|l . Substituting equation Ijll^ into equation Q, we obtain 

{Cy^u)„)ab dx" A dx^ = d[u^^{{l - e„)Qn + e„flp}] A dD,n = , (12) 

{CvpOJp)ab dx" A dx'' = d[ti7^{(l — ep)flp + epD.„}] A dQ.p = , (13) 

where zu is the distance from the rotation axis, defined hy zu — r cos 9. Thus, these vorticity conservation equations are 
automatically satisfied if the rotation laws for two superfluids are given as follows: 

- £„)il„ + Eniip} = j„{^„) , - £p)Qp + EpQn} = jp(^p) , (14) 

where j„{^l„) and jp(r2p) are arbitrary functions of Sl„ and ilp, respectively. In other words, the rotation laws of superfiuid 
stars must be strongly restricted because the rotation velocities Sl„ and i}p may depend on the entrainment functions £„ and 
Bp, whose functional forms are in general determined with the matter distributions of stars. When the effect of the entrainment 
between two superfiuids does not operate, that is, the case of a = or of S7„ = Sip, the vorticity conservation equations become 
relatively simple as follows 

dro A drj„ = , d'cuAdQ,p = 0. (15) 

The functions S7„ and flp must not therefore depend on 2 = r sin 61 in this situation. In other words, we can choose the rotation 
velocities n„ and Qp freely as long as Q„ — Q„{-ou) and Qp = i}p{-uj) are fulfilled. Note that this condition is the same as that of 
barotropic ordinary fluid stars. For simplicity, in this paper, we assume that both the neutrons and the protons are uniformly 
rotating, i.e. Q„ = const and Qp — const, because vorticity conservations are automatically satisfied for uniformly rotating 
configurations. Assuming the uniform rotations and integrating the Euler equations Q and ((TJ, we can obtain integrated 
equations for hydrostatic equilibrium states, given by 
2 2 

fiu + $ - —ni = Cn, fip + $-—np = Cp, (i6) 

where C„ and Cp are integral constants. Here, we have used the relationships 
Cv„{Vn + £nA)adx'' = UJ^ {{I — e„)Q„ + e„Qp} dD,n = , 

Cvp{vp ~ £pA)adx° = ~ ep)^lp + £pi^„} dQp ^ . (17) 

We are interested in configurations in which two superfluids are rotating rapidly with different angular velocities, f2„ 7^ Sip. 
We do not therefore make use of any slow-rotation approximation such as the Chandrasekhar-Milne expansion. We however 
assume the difference between rotation velocities of two superfiuids to be very small. We can then treat the effect of the 
angular velocity difference Qp — fl„ on the stellar structure as a perturbation to rapidly rotating stars whose rotation rates 
of two superfiuids are exactly the same. For superfiuid neutron stars, although the angular velocity S7„ may differ from the 
angular velocity flp because the interaction between two superfluids may be very weak due to the superfluidity, it is believed 
that this difference between two angular velocities is small in comparison with the averaged angular velocity of the star. Thus, 
this treatment could be appropriate for neutron star models with superfluidity. 

We expand physical quantities appeared in equations JHJ and llbl l up to the first order in ($!„ — f2p)/(|f2„| + \Qp\) as 
follows: 

n„ ^ n{i + SQr,) , np^n{i + 5Qp), (is) 

Pn = PnO + Sp„ , Pp = PpO + 5pp , (19) 
Pn = AnO + Sp„ , Mp = ApO + Spp , (20) 

<l> = <l>o + 5$, C„ = C„o + <5C„ , Cp^Cpo + 5Cp, (21) 
where jlno, Ppo, Sjln, and 5 ftp are given in terms of partial derivatives of the internal energy density £, Spn, and Spp by 
d£ . d£ 

PnO = -5 , PpO = , (22) 

opn app 



d^£ d^£ d'^£ d^£ 

dpi apndpp dpi oppOpn 
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Here, quantities following "5" denote perturbed ones. Note that since we consider the accuracy up to the first order in the 
rotation velocity difference, no entrainment effects appear in equations we solve. Substituting those expanded quantities into 
equations ((HJ and If 6II . we obtain the zero-th order equations, given by 

fiuo + $0 - — = C„o , (24) 

2 

ApO + 'l'o-^f^'^Cpo, (25) 

V„V"$o = 47rG(p„o + Ppo) . (26) 
Equations 1241 and I25I I lead to a relationship between chemical potentials of the neutron and the proton 

jj.nO — jjpQ = CnO ~ CpO . (27) 

In this paper, we take CnO and Cpo to be Cno = Cpo, assuming that the two superfluids are in chemical equilibrium in the 
unperturbed state. Thus, our basic equations for obtaining unperturbed rapidly rotating stars can be written by 

^2 

/inO + "I'd ~ CnO , jipO ~ flnO ■ (28) 

On the other hand, the equations for the perturbed quantities are given by 

5p,n + 5* - Vj'^n'^SQn = SCn , (29) 

5fj,p + 5$- vj'^n'^snp = sCp , (30) 

VaVS^^AnGiSpn + Spp). (3f) 

In order to solve the Poisson equations II26|I and l|8f |l . boundary conditions are required. The appropriate boundary conditions 
at spatial infinity and at the center of the star are given by 

$0^0, (5$^0, as r^oo, 9^$o = , drS^ = , at r = . (32) 
2.3 Equation of state 

In order to obtain equilibrium configurations, the internal energy density £ must be specified. In this paper, we employ the 
same analytical internal energy density as that used in Prix et al. (2002a), which is given by 

£ = 2kxp{l - Icpia + 1)} ^^"^^ + {1 - + '^(1 - 2a;p)}Pp - 2aXpPnpp] , (33) 

where Xp, a, and k axe constants. Note that the entrainment that is included in the equation of state of Prix et al. (2002a) 
is neglected because the entrainment appears in second order in 5Qn or SQp. This internal energy density is a natural 
generalization of an A'^ = f polytropic equation of state for a barotropic ordinary fluid into the two superfluids because it is 
written in the general quadratic form of Pn and Pp, given by 

£ ~ t^nnPn ^" '^f^npPnPp + ^ppPp 5 (^4) 

where Knp, and npp are constants. The chemical potentials for the internal energy density I33II are then given by 

ftn = 77- — —TyriPn-app), 

k{l - Xp[a + 1)1 

ftp = ^^^^^ _ ^^^^ ^ [{f -Xp + a{l - 2xp)}pp - axppn] . (35) 
On the other hand, p„ and pp can be written in terms of /i„ and pp as 

Pn = ^^[{1 - Xp + a{l - 2Xp)}pn + ^Xpflp] , Pp = [apn + Pp) , (36) 

O" + i O" + i 

which lead to 

p = Pn + Pp = k{{l - Xp)p„ + XpPp} . (37) 

Let us consider the situation where the neutrons and the protons are in chemical equilibrium, which is given by the condition 
p.n ~ Pp- From equation H35^ . we can obtain 

XpPn = (1 — Xp)pp . (38) 
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This equation means that Pn/pp = constant inside the star and that Xp represents the proton fraction when chemical 
equilibrium between the neutrons and the protons is achieved. On the other hand, the parameter cr of the internal energy 
density 13311 is interpreted as the so-called "symmetry energy" term (Prix et al. 2002a; Prakash, Lattimer, & Ainsworth 1988). 
The appropriate range for a might be cr £ ( — 1, 1) and we consider only three values for a, i.e. a = —0.5, 0, 0.5, in this paper 
(Prix et al. 2002a). 



2.4 Unperturbed state: rapidly rotating stars 

For numerical calculations, it is convenient to introduce non-dimensional physical quantities as follows: 

r — rof , pn — Papn , Pp — pOpp , Pn = POPn , Pp = POPp , 

n„ = ^/I^vGp^Cln , ilp = y^InGfX) Ctp , $ = 47rGpo'"o'i> , k = kok , (39) 
where quantities with hat are non-dimensional ones, and ro and are defined by 



ro = y/l/i^nGko) , po = Po/ko , (40) 

where k is determined so as to be fmax = 1 for the unperturbed star, and fmax is the largest distance from the stellar surface 
to the center of the unperturbed star. Here, ro and po can be given freely because of the polytropic equation of state. 

Since, for unperturbed states, chemical equilibrium between two superfluids are assumed, the chemical potential pno can 
be written as pno ~ k~^{l ~ Xp)~^pno ~ k~^po by virtue of the analytical equation of state I35L where po is the total mass 
density. Then the master equations for unperturbed stars are reduced to 

k-^po + $0 - ^ = CnO , (41) 

1 f Mn^^s^,^ 

r^dr'hn{r,r) / sm 9' dO' P2„{cos 6') po{f ,0') , (42) 



47r / r - r' 



OC „c 

^ Pan (cose) / 

n=0 "'0 







where 



I f' 



1 ^{frY" for r'>f, 

and P2n (cos 9) are the Legendre polynomials. Here, the gravitational potential $o has been written in the integral representa- 
tion, in which a proper Green function is employed to include physically appropriate boundary conditions both at f = and at 
spatial infinity 1321 . In this integral representation, equatorial symmetry of the matter distribution has been assumed because 
we are concerned about stars having this equatorial symmetry. Note that these equations are exactly the same as those of 
polytropic ordinary fluid stars with polytrope index A'^ = 1. Thus, solutions to equations H41|l and l|42^ can be obtained by the 
so-called Hachisu's Self Consistent Field scheme (HSCF), once the axis ratio, rp, is given (Hachisu 1986). Here, rp is defined 
by 

^p — rmin/rniax — rniin , (^^) 

where fmin is the minimum distance from the stellar surface to the center of the star. Note that solutions po, $o, and C„o are 
universal functions in the sense that those are independent of the proton fraction Xp. After getting a solution po, and giving 
the proton fraction Xp, we can obtain the mass densities and chemical potentials of two superfluids through the formulas 

PnO = (1 ^ ^p)PO , PpO = ^pPO , PnO = PpO = Pok~^ ■ (45) 

The surface of the star R{9) for unperturbed stars is defined by 
P{R{9),9)=0, (46) 
where P is the generalized pressure, which is given, for our internal energy density, by 

P = — —[^ppl -I- {1 - a;p + o-(l - 2xp)}pl ~ 2axpp„pp] . (47) 

2fca:p{l - Xp(a + 1)| 

Pressure Po for the unperturbed state can be explicitly written as 

n POPO .2 POPO .2 

'^'=2H1-XpY'-' = ^''- ^''^ 
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The non-dimensional stellar surface Ro{6) for an unperturbed star is then given by a solution of the algebraic equation 
2.5 Effect of the rotation velocity difference on the stellar structure 

In this paper, we consider perturbed states whose densities of two superfluids p„ and pp at the center of the star have the same 
values as those of unperturbed states. Because of our equations of state HH^fl . this requirement is equivalent to the assumptions 
of Sjin = = at the center of the star. Thus, the integral constants in equations 1)29^ and may be determined so as 
to be 

5C„ = SCp = <5$(r = 0, 6*) . (49) 

For the analytical internal energy density . perturbations of the total mass density 5p is given in terms of perturbed 
chemical potentials 5fi„ and 5 ftp as 

5p — 5pn + 5pp = k {(1 — Xp)5jin + XpSfip} . (50) 
From equations 129^ and H3()|l . thus, we obtain 

k'^Sp + (5$ - m^Cl^{{l - Xp)5fl„ + XpSflp} = SCn , (51) 

where the relations of the integral constants l|49|l have been assumed. Note that the perturbed total density is dependent on 
the proton fraction Xp but not on the symmetry energy parameter a. From equation l|51|l . it is found that perturbations 5p, 
5$, and 5Cn can be represented in terms of three functions independent of SQ„, 5Qp, and Xp as follows: 

Sp = {{1- xp)5nn + xp5np}Sp, s$ ^ {{1 - Xp)5nn + xpSQpjs^ , 
5c„ = {{1 ~ xp)5nu + xpSnpjsCn , (52) 

where the three functions 5p, 5^, and 5Cn are the solutions of the equations 

k'^Sp + S^ -m^Cf ^ 5Cn, (53) 

VaV^S^^Sp. (54) 

Note that three quantities 5p, 5$, and 5Cn are universal functions in the sense that those are independent of the parameters 
Xp and a, and that those are dependent only on the structure of the unperturbed star. With equations 12911 and 130L the 
chemical potentials of two superfluids are given by 

5p„ = {xp{S^ ~ SCn) + k~^5p}5nn ~ Xp{S^ ~ SC„)SQp , 

Spp = {l-Xp){SCn-S^)Snn + {{l-Xp){5^-5Cu)+k-^5p}Snp. (55) 
Perturbations of the mass densities, on the other hand, can be written in terms of Spn and Sfip as 

Sp„ = — ^— [{1 - Xp + a(l - 2xp)}Sfl„ + aXpSfip] , 5pp = (^aSfin + Sfip) , (56) 

(T + i (7 + i 

Note that, for our internal energy density, the chemical potentials do not depend on the symmetry energy parameter a, while 
the mass densities do. The perturbations of the generalized pressure SP is given by 

5P = popok^^ po{{l - Xp)Sp„ + XpSpp} . (57) 
If we write the stellar surface with accuracy up to the first order of the rotation velocity difference as 

R{e) ^ Roie) + 5R{e) , (58) 
the first order corrections for the stellar surface 5R{9) can be due to equation 14611 written in terms of Sfi„ and Sfip as 
SR{e) = (1 - Xp)5R„ + XpSRp , (59) 
where SRn and SRn are defined by 

5R,. = -^{f = Ro{e),6), 5Rp = -^{f = Ro{e),6). (60) 

OrpO OrpO 

Because fin{Ro + SRn) = and fip{Ro + SRp) = are satisfied, two surfaces of r = _Ro + SRn and f = Ro + SRp can 
be interpreted as the surfaces of the neutron and proton superfiuids, respectively. Although we can define the respective 
fiuid surfaces of two superfiuids as the zero-density surfaces (Prix 1999; Prix et al. 2002a), we make use of the definition 
with the chemical potential for the surfaces because the definition with the chemical potentials is considered to be natural 
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generalization of the surface definition for ordinary fluid stars, for which the surface is defined as a zero-pressure surface not 
as a zero-density surface. Note that as we can see from equation (35), in general, "zero-chemical potential surfaces" do not 
coincide with "zero-density surfaces". In other words, the densities do not necessarily vanish on the zero-chemical potential 
surfaces because the equi-potential surfaces in general incline to the equi-density surfaces in the two-fluid model. 

In our numerical procedure of solving equations 1)53^ and 1)54^ . in order for the Poisson equation H54|l to satisfy the 
boundary condition H32|l explicitly, equation l|54|l is converted into the integral representation, given by 

5$ = -^P2n(cos6») / sine'de' r''drf2nif,r)P2n{co8e')Sp{r,e'). (61) 

In this paper, equations 15311 and 16111 are numerically solved with a variant of the so-called Self Consistent Field scheme 
(Ostriker & Mark 1968). To obtain solutions, we follow the following steps: i) By assuming an initial guess for 5p, compute the 
gravitational potential through the two-dimensional integration H61|l . ii) Determine the value of parameter 5C„ from equation 
igg. iii) By using the obtained 5$ and SCn, solve equation l|53|l for the perturbed density Sp. iv) Compare the obtained 
perturbed density with the one used in obtaining the perturbed potential. If the relative changes of Sp are less than 10"** at 
all grid points, then the obtained perturbed density distribution is considered to be a converged solution. If the condition for 
the relative changes is not satisfied, go back to step ii) and there the obtained perturbed density is treated as a new initial 
guess for Sp for the next iteration cycle. 

In actual computations, we make use of equidistantly spaced discrete meshes in the radial direction (0 < f < fmax = 1). 
In order to calculate integrations in r direction, we employ a classical trapezoidal rule. As for the angular variable 9, we take 
the angular grid points at pi = cos^i, where pi's are zeros of the Legendre polynomial of order 2L — 1, i.e. P2L-i(pi) = 0, 
and there are L grid points for Q < 9 < ti/2. We employ the Caussian quadrature (see, e.g., Abramowitz & Stegun 1964) to 
evaluate integrations in the 9 direction. In the present investigation, we take the number of mesh points to be 500 x 25 (r x9). 
The Legendre polynomials in equations l|42|l and H61^ are added up to Pie (p) . 



3 NUMERICAL RESULTS 

In this paper, we assume the proton fraction to be Xp — 0.1 because a typical proton fraction in the cores of old neutron stars 
is expected to be around Xp — 0.1. For the parameter of the symmetry energy, a, since we have little information about the 
range of a, we investigate three cases, i.e. a = —0.5 , , 0.5, in order to examine the effect of the symmetry energy on the 
stellar structures. Note that the total density of superfluids and the chemical potentials are dependent on the proton fraction 
but not on the symmetry energy. We exhibit the results for two cases of {Si^n,Si^p) = (1,0) and (5Q,ri,S^lp) — (0, 1) because 
solutions for other values of {SQn, 5Q,p) are expressed in terms of linear superpositions of those two solutions. 

First, let us consider the unperturbed stars. In our unperturbed states, two superfluids are assumed to be in the same 
rotational motion and in chemical equilibrium. As mentioned in the last section, the structures of unperturbed stars are, 
therefore, almost the same as those of uniformly rotating = 1 polytropic stars. Thus, there are no new features for the 
unperturbed stars. For reader's convenience, however, we exhibit some results for the unperturbed stars. In Figures 1 through 
4, we show the non-dimensional fundamental quantities of the unperturbed stars, the axis ratios, Vp, the total masses, M, 
the total moments of inertia, /, and the ratios of the rotational energy to the absolute value of the gravitational energy, 
r/|Vy|, as functions of the angular velocity, Q, where Q is the non-dimensional angular velocity of the star, defined by 
n = Q,/ {GM / P?Y^^ , and R denotes the stellar radius on the equatorial plane. Here, the total mass, M is defined by 

M = j pod'r = yPo^oM = MoM . (62) 

The total moments of inertia, I, is given by 

Pouj^d^r = 47rporo/ = IqI . (63) 

The rotational energy, T, and the gravitational energy, W, are, respectively, defined by 

T =\ I povj^^'d'v, W = [ po'i'od^r . (64) 



27 2 

Since Pn/pp = const for the unperturbed stars, we can write the masses, M„, Mp, and the moments of inertia, I„, Ip, for the 
neutrons and the protons as 

M„ = (l-Xp)M, Mp=XpM, In^{l-Xp)I, Ip = XpI. (65) 

In order to check our numerical code, we have calculated a solution for a slowly rotating star with the axis ratio Vp = 0.992, 
whose angular velocity is H = 0.10298, and compared them with the analytical solution obtained by Prix et al. (2002a), in 
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which the structure of a slowly rotating superfluid star was examined analytically with the slow rotation approximation. The 
perturbed densities obtained from two methods are shown as functions of r in Figure 5. The solid lines show the analytical result 
for 5pn with the slow rotation approximation, and the solid squares our numerical result without assuming the slow rotation 
approximation. The parameters of the model shown in Figure 5 have been chosen to be a = —0.5 and {5Qn, Silp) = (1, 0). As 
seen from this figure, our result is in good agreement with the analytical result of Prix et al. (2002a) as long as the angular 
velocity of the star is small enough. 

In Figures 6 through 11, typical density distributions for the neutrons and the protons, Sp„ and Spp are represented 
for three values of the symmetry energy parameter, a = —0.5, 0, 0.5. For the models shown in Figures 6 through 11, the 
axis ratio rp for the unperturbed star has been taken to be rp = 0.668, and the corresponding angular velocity is given 
by n = 0.8037. Figures 6 through 8 show the density distributions for {SQn, SQp) = (1,0), and Figures 9 through 11 for 
{Sfln, 5Qp) — (0, 1). In those figures, the perturbed densities are shown versus the radial coordinate r for six different values of 
9i, in which the longest curve corresponds to the result on the equator, 9 = 7r/2, and di's decrease for each successively shorter 
curve toward the value = on the symmetry axis. Comparing Figure 6 (left panel) with Figure 5, in which distributions 
of Sp„ are shown for the same parameter as those in Figure 6 except for the axis ratio Vp or for the rotation velocity H, we 
observe that the basic qualitative properties of the density perturbations for two superfluids do not depend on the value of Q 
so much. We do not therefore show the density distributions for other rotation rates Cl in this paper. From Figures 6 through 
11, we can see that the amplitudes of Sp„ are much larger than those of Spp for the solutions of {SQ„, S^lp) — (1, 0), while the 
amplitudes of Sp„ is smaller than those of Spp for the case of (SQn, Sflp) = (0, 1). As shown in Figures 9 through 11, however, 
for the solutions of {SQn, SQp) — (0, 1), difference of the amplitude between Spn and Spp is not so large. This is because the 
proton fraction is not very large in the interior. Similar behaviors were found by Prix et al. (2002a). It is also found from 
Figures 6 through 11 that Spp is strongly dependent on the values of a for the {SQn, SQp) = (1,0) case, while Spn depends 
strongly on a for the {SQ„,5np) = (0, 1) case. This means that the significant parameters to model a superfluid neutron star 
are appropriate combinations of the symmetry energy parameter a and the velocity difference between two superfluids. 

In Figures 12 and 13, the perturbed axis ratios are plotted as functions of the angular velocity f2. Here, the functions, 
Svp, Svp^n, and Svp^p can be interpreted as the perturbations of the axis ratios of the whole star, the neutron superfluid, and 
the proton superfluid, respectively, and are defined by 

Srp = SR{e = 0)-rpSR{e^Tr/2), Srp^„ ^ SR„{9 ^ 0) - rpSRn{9 = -k/2) , 
5rp,p = 5Rp{9 = 0)-rpSRp{9 = -K/2). (66) 



The results for {SQn, Sflp) = (1,0) and {S^ln, SQp) — (0,1) are shown in Figures 12 and 13, respectively. Note that in these 
figures, the perturbed axis ratios are shown for < H < 0.9 because they diverge as Q goes to its maximum value. These figures 
illustrate how the surface of the star is deform due to the effect of the rotation velocity difference between two superfluids. 
It is found that the change of the axis ratio of the rotating component with the slower angular velocity becomes positive 
when the angular velocity Q, becomes larger than some critical value. In other words, the rotating component with the slower 
angular velocity tends to become prolate when a star rotates very rapidly. 

In Figure 14, the perturbed mass is shown as a function of H. Here, the perturbed total mass SAl is deflned as 

SM = — ^ / Spd^r . (67) 
4ttM J 

The solid curve and the dashed curve show SM for the solutions with {SQn, SQp) = (1, 0) and {SCln, SQp) = (0, 1), respectively. 
In this paper, we consider the perturbed states whose central density is the same as that of the unperturbed states. Thus, the 
mass of the perturbed stars can change. In Figures 15 and 16, we plot the perturbed moments of inertia for the neutron and 
the proton as functions of the angular velocity Q. Here, the perturbed moments of inertia for the neutron and the proton, SI„ 
and Sip, are defined as 

SI„ — ^ / Spn'co^ d^Y , Sip — i / SppUp' (^v . (68) 

^■Krl{\ - Xp)I J A-Kr^XpI J 

Figures 15 and 16 show the perturbations of the moments of inertia for {SflnjSflp) = (1, 0) and {SQ„, Sflp) = (0, 1), respectively. 
In these figures, results for three different values of a, a = —0.5, 0, 0.5 are displayed. Figures 15 and 16 show that the perturbed 
moments of inertia for the proton are strongly dependent on the symmetry energy parameter a for {5Qn,S^lp) — (1,0) case, 
while the perturbed moments of inertia for the neutrons strongly depend on a for {Sfl„, SQp) — (0, 1) case. This is consistent 
with that of the properties of the density distributions for the neutrons and the protons. It is noted that amplitudes of SI„ 
and Sip become large as the angular velocity Q is increased because the perturbations are roughly proportional to Ct^ due to 
equation 1531 . 
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4 SUMMARY AND DISCUSSION 

In this paper we have developed a formulation for constructing and examining rapidly rotating neutron stars that contain 
two superfluids, taking account of the effect of the rotation velocity difference between two superfluids. We assumed neutron 
stars to be composed of the superfluid neutrons and the mixture of the superfluid protons and the normal fluid electrons. 
To describe Newtonian dynamics of the two superfluids, the Newtonian version of the two-fluid formalism developed by Prix 
(2002) was employed. In this paper, we considered the situation where two superfluids rapidly rotate around the same rotation 
axis, but the rotation velocity difference between two superfluids is very small. We then treated the effect of the rotation 
velocity difference on the equilibrium configurations as a small perturbation to rapidly rotating superfluid stars whose rotation 
velocities of two superfluids are the same. We derived basic equations for perturbations on the structure of rapidly rotating 
superfluid stars due to the rotation velocity difference between two superfluids. Assuming the superfluids to obey a simple 
analytical equation of state used by Prix et al. (2002a), we obtained numerical solutions for the perturbed quantities and 
found that the density distributions of the superfluids are strongly dependent of the symmetry energy parameter a, which 
appears in the analytical equation of state. Similar properties were found in Prix et al. (2002a). It was also found that if 
Prix et al. (2002a) 's analytical equation of state is assumed, the perturbations can be represented in terms of the universal 
functions that are independent of the parameters of the equation of state. 

Although we only considered one special equation of state for superfluids in the present investigation, the formalism we 
derived in this paper is straightforwardly applicable to general equations of state. We treated the superfluids inside neutron 
stars in the framework of Newtonian dynamics. The effect of general relativity must not be however neglected for the structures 
of neutron stars because neutron stars are quite compact in a sense that a general relativistic effect can be expressed by the 
factor GM/(?R whose typical value is ~ 0.2 for neutron stars, where c is the speed of light. It is straightforward to extend 
the present formulation to general relativistic configurations and we will do it in the future. Recently, Prix et al. (2002b) have 
obtained rapidly rotating relativistic superfluid stars whose rotation velocity of the neutrons differs from that of the protons. 
They did not assume that the rotation velocity difference between two fluids is very small. In other words, their method 
can be applied to equilibrium configurations with any rotation velocities. Yet, our perturbation method developed in this 
investigation has an advantage in studying structures of real neutron stars because it would give results with higher accuracy 
by taking account of the smallness of the velocity difference directly, as long as the rotation velocity difference between two 
fluids in real neutron stars is very small. A solid crust is believed to exist near the surface of cold and old neutron stars, 
and can have a significant influence on the structure of the equilibrium states if the solid crust is not in a strain-free state in 
the equilibrium states of rotating neutron stars. Therefore, an investigation of the effect of the solid crust on the equilibrium 
states remains as a challenging problem in the future too. 

Acknowledgments: S.Y. acknowledges financial support from Fundagao para a Ciencia e a Tecnologia (FCT) through project 
SAPIENS 36280/99. 
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Figure 1. Axis ratio rp for the unperturbed stars is given as a function of (i 
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Figure 2. Same as Figure 1 but for the mass M 
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Figure 4. Same as Figure 1 but for the ratios of the rotational energy to the absolute value of the gravitational energy T/|W| 
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Figure 6. Perturbed densities for the neutron (left) and the proton (right) are given as functions of r for six different values of 6. The 
model parameters are Q = 0.8032, a = —0.5, and {5Qn,SClp) = (1, 0). 
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Figure 10. Same as Figure 9 but for u = 0. 
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Figure 11. Same as Figure 10 but for a = 0.5. 
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Figure 12. Perturbed axis ratios Srp, Srp^n, and Srp^p are given as functions of Q. The solid, dashed, and dotted curves denote Sr^ 
Srp^n, and Srp^p, respectively. The model parameter is {SQn,Siip) = (IjO). 
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Figure 13. Same as Figure 12 but for {Sdn,SClp) = (0, 1). 
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Figure 14. Perturbed mass 5M is given as a function of Q. The solid and dashed curves denote results for {SQn,SQp) = (1,0) and 
{SQn,SOp) = (0, 1), respectively. 
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Figure 15. Perturbed moments of inertia for the neutron and tlie proton iln, 6Ip are given as functions of O. The model parameter is 
(5n„,<5np) = (1,0). 
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Figure 16. Same as Figure 15 but for ((5f2„,(5np) = (0,1). 



